﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "gaus_int")]
    public static unsafe double gaus_int(int n, IntPtr js_ptr, IntPtr s_xa_ya_mn_ptr, IntPtr f_xa_n_ptr)
    {
        int* js = (int*)js_ptr.ToPointer();
        s_xa_ya_mn = Marshal.GetDelegateForFunctionPointer<delegatefunc_xa_ya_mn>(s_xa_ya_mn_ptr);
        f_xa_n = Marshal.GetDelegateForFunctionPointer<delegatefunc_xa_n>(f_xa_n_ptr);

        return gaus_int(n, js);
    }

    /// <summary>
    /// 计算多重积分的高斯法
    /// s(double [],double [],int,int)计算各层积分上下限（要求所有的上限>下限）的函数名。
    /// f(double[], int)计算被积函数值f(X)的函数名。
    /// double gaus_int(int n, int js[], void (*s)(), double (*f)(int,double []))
    /// </summary>
    /// <param name="n">积分重数。</param>
    /// <param name="js">js[n] js[k]表示第k层积分区间所划分的子区间个数。</param>
    /// <returns>函数返回积分值。</returns>
    public static unsafe double gaus_int(int n, int* js)
    {
        int m, j, k, q, l;
        double p;
        double[] t = new double[5] {
            -0.9061798459,-0.5384693101,0.0, 0.5384693101,0.9061798459
        };
        double[] c = new double[5] {
            0.2369268851,0.4786286705,0.5688888889,0.4786286705,0.2369268851
        };
        int* si = stackalloc int[2 * (n + 1)];
        double* x = stackalloc double[n];
        double* a = stackalloc double[2 * (n + 1)];
        double* b = stackalloc double[n + 1];
        double* y = stackalloc double[2];

        m = 1;
        l = 1;
        a[n] = 1.0;
        a[2 * n + 1] = 1.0;
        while (l == 1)
        {
            for (j = m; j <= n; j++)
            {
                // 计算j-1层积分区间的上下限y[1]与y[0]
                s_xa_ya_mn(x, y, j - 1, n);
                a[j - 1] = 0.5 * (y[1] - y[0]) / js[j - 1];
                b[j - 1] = a[j - 1] + y[0];
                // 高斯点
                x[j - 1] = a[j - 1] * t[0] + b[j - 1];
                a[n + j] = 0.0;
                // 这是j-1层积分的第1个子区间的第1个高斯点
                si[n + j] = 1;
                si[j - 1] = 1;
            }
            // 从最内层积分开始
            j = n;
            q = 1;
            while (q == 1)
            {
                //取j-1层积分区间当前子区间上的高斯点序号
                k = si[j - 1];
                //计算高斯点上的被积函数值
                if (j == n) p = f_xa_n(x, n);
                else p = 1.0;
                a[n + j] = a[n + j + 1] * a[j] * p * c[k - 1] + a[n + j];
                //置j-1层积分区间当前子区间的下一个高斯点序号
                si[j - 1] = si[j - 1] + 1;
                //j-1层积分区间当前子区间上的高斯点全部计算完
                if (si[j - 1] > 5)
                {
                    //j-1层积分区间的所有子区间考虑完
                    if (si[n + j] >= js[j - 1])
                    {
                        // 考虑前一层的积分区间
                        j = j - 1;
                        q = 1;
                        // 已到最外层
                        if (j == 0)
                        {
                            p = a[n + 1] * a[0];
                            return (p);
                        }
                    }
                    else
                    {
                        // j-1层积分区间还有子区间
                        // 置j-1层积分区间的下一个子区间
                        si[n + j] = si[n + j] + 1;
                        b[j - 1] = b[j - 1] + a[j - 1] * 2.0;
                        // 这是j-1层积分区间当前子区间的第1个高斯点
                        si[j - 1] = 1; k = si[j - 1];
                        x[j - 1] = a[j - 1] * t[k - 1] + b[j - 1];
                        // 这是最内层
                        if (j == n) q = 1;
                        // 这不是最内层
                        else q = 0;
                    }
                }
                else
                {
                    // 计算j-1层积分区间当前子区间上的下一个高斯点
                    k = si[j - 1];
                    x[j - 1] = a[j - 1] * t[k - 1] + b[j - 1];
                    // 这是最内层
                    if (j == n) q = 1;
                    // 这不是最内层
                    else q = 0;
                }
            }
            m = j + 1;
        }
        return (0.0);
    }

    /*
    // 计算多重积分的高斯法例
      int main()
      { 
          int js[3]={4,4,4};
          void gauss(int,int,double [],double []);
          double s,gausf(int,double []);
          s=gaus_int(3,js,gauss,gausf);
          cout <<"s = " <<s <<endl;
          return 0;
      }
    // 计算各层积分上下限
      void gauss(int j, int n, double x[], double y[2])
      { 
          double q;
          n=n;
          switch (j)
          { 
            case 0: { y[0]=0.0; y[1]=1.0; break;}
            case 1: { y[0]=0.0; y[1]=sqrt(1.0-x[0]*x[0]); break;}
            case 2: { q=x[0]*x[0]+x[1]*x[1]; y[0]=sqrt(q);
                      y[1]=sqrt(2.0-q); break;
                    }
            default: { }
          }
          return;
      }
    // 计算被积函数值
      double gausf(int n, double x[])
      { 
          double z;
          n=n;
          z=x[2]*x[2];
          return(z);
      }
    */
    /*
    // 计算多重积分的高斯法例
      int main()
      { 
          int js[3]={10,10,10};
          void gauss(int,int,double [],double []);
          double s,gausf(int,double []);
          s=gaus_int(3,js,gauss,gausf);
          cout <<"s = " <<s <<endl;
          return 0;
      }
    // 计算各层积分上下限
      void gauss(int j, int n, double x[], double y[2])
      { 
          n=n;
          switch (j)
          { 
            case 0: { y[0]=1.0; y[1]=2.0; break;}
            case 1: { y[0]=1.0; y[1]=2.0; break;}
            case 2: { y[0]=1.0; y[1]=2.0; break;
                    }
            default: { }
          }
          return;
      }
    // 计算被积函数值
      double gausf(int n, double x[])
      { 
          double z;
          n=n;
          z=x[0]*x[0]+x[1]*x[1]+x[2]*x[2];
          return(z);
      }
    */
}

